Edge-preserving enhancement of seismic images by nonlinear anisotropic diffusion

ABSTRACT

rocessing a seismic image comprising the steps of(a) obtaining a two or three dimensional initial image data set, wherein each element u(m=0) of the data set is the initial image intensity of a point of the image;(b) calculating for each point the partial derivatives of u(m) in n directions to obtain n derivative data sets of the partial derivatives u xi ; (c) calculating for each point a symmetric n×n structural matrix S, wherein the elements s pq  equal u xp ·u xq , and an n×n diffusion matrix D, wherein the elements d ij  are a function of the elements s pq ; (d) calculating for each point u(m+1) from u(m) using the following equation u(m+1)=u(m)+c.div(ε.Dgrad(u(m))), wherein c is a predetermined constant, 0≦c≦1, and wherein ε is a scalar, 0≦ε≦, wherein ε is close to zero when near an edge and close to 1 when far away from an edge ( 10, 11, 12, 14, 15 ); and(e) repeating steps (b) through (e) M times to obtain the processed image.

[0001] The present invention relates to a method of processing an image,for example a seismic image, applying a diffusion process. Such methodsare known, see for example the thesis ‘Anisotropic diffusion in imageprocessing’ by J Weickert, January 1996, pages 41-44.

[0002] Known is a diffusion method that comprises the steps of

[0003] (a) obtaining an n-dimensional initial image data set, whereineach element u(m=0) of the data set is the initial image intensity of apoint of the image, and wherein n=2 or n=3;

[0004] (b) calculating for each point the partial derivatives of u(m) inn directions to obtain n derivative data sets of the partial derivativesu_(xi);

[0005] (c) calculating for each point a symmetric n×n structural matrixS, wherein the elements s_(pq) equal u_(xp)·u_(xq);

[0006] (d) calculating for each point an n×n diffusion matrix D, whereinthe elements d_(ij) are a function of the elements s_(pq);

[0007] (e) calculating for each point u(m+1) from u(m) using thefollowing equation u(m+1)=u(m)+div(Dgrad(u(m))); and

[0008] (f) repeating steps (b) through (e) M times to obtain theprocessed image.

[0009] Unless otherwise specified, the following holds throughout thespecification and the claims:

[0010] i, j, p, q=1, . . . , n ;

[0011] u_(xi)=∂u/∂x_(i).

[0012] An element of the image data set is also referred to as a pointof the image.

[0013] In the known method an initial image is processed to obtainsuccessive images, wherein each image is obtained from the previous oneby a diffusion process. The diffusion method can as well be described bythe${{diffusion}\quad {equation}\quad \frac{\partial u}{\partial t}} = {{{div}\left( {D \cdot {\nabla\quad u}} \right)}.}$

[0014] There are four types of diffusion methods. At first there is theisotropic diffusion in which the diffusivity D is a scalar and theanisotropic diffusion in which D is a matrix. Then for each of these twotypes, the diffusion can be linear, in which D is not a function of u orthe diffusion can be non linear in which D is a function of u.

[0015] This method can as well be applied to processing a seismic image,in that case the image intensity is the amplitude of the seismic signal.

[0016] An advantage of the diffusion method is that noise is suppressed.A further advantage, which is particularly relevant to processing aseismic image is that after a few cycles the geological structure ishighlighted more clearly. However, a problem with this method ispreserving edges.

[0017] Edge preservation is of great interest in the interpretation ofseismic images, which play an important role in the study of undergroundformations, and in particular in the study of underground formationsthat contain hydrocarbons. Moreover there is a great interest in edgepreserving techniques, which allow highlighting of faults while removingnoise from the seismic data.

[0018] For edge preservation in combination with isotropic diffusionprocessing is proposed a method based on the Perona-Malik model (see thethesis of Weickert, page 10) of which the diffusion equation is${\frac{\partial u}{\partial t} = {{div}\left( {{g\left( {{\nabla\quad u}}^{2} \right)}{\nabla\quad u}} \right)}},$

[0019] wherein the non-linear diffusivity is g(|∇u|²)=(1+|∇u|²/λ²)⁻¹,the scalar λ being a predetermined edge preservation parameter.

[0020] Applicant proposes a simple addition to the known method that hasproved to be an effective way of preserving edges while capable ofremoving noise. Moreover Applicant proposes such a method that isapplicable to anisotropic diffusion.

[0021] To this end the method of processing an image according to thepresent invention comprises the steps of

[0022] (a) obtaining an n-dimensional initial image data set, whereineach element u(m=0) of the data set is the initial image intensity of apoint of the image, and wherein n=2 or n=3;

[0023] (b) calculating for each point the partial derivatives of u(m) inn directions to obtain n derivative data sets of the partial derivativesu_(xi);

[0024] (c) calculating for each point a symmetric n×n structural matrixS, wherein the elements s_(pq) equal u_(xp)·u_(xq);

[0025] (d) calculating for each point an n×n diffusion matrix D, whereinthe elements d_(ij) are a function of the elements s_(pq);

[0026] (e) calculating for each point u(m+1) from u(m) using thefollowing equation u(m+1)=u(m)+c.div(ε.Dgrad(u(m))), wherein c is apredetermined constant, 0≦c≦1, and wherein ε is a scalar, 0≦ε≦1, whereins is close to zero when near an edge and close to 1 when far away froman edge; and

[0027] (f) repeating steps (b) through (e) M times to obtain theprocessed image.

[0028] In this way edge preservation is introduced in anisotropic,non-linear diffusion.

[0029] Suitably step (c) comprises calculating for each point asymmetric n×n structural precursor matrix S⁰, wherein the elements s⁰_(pq) equal u_(xp)·u_(xq); creating (½)n(n+1) data sets, wherein theelements of the first data set are s⁰ ₁₁ pertaining to each point, theelements of the second data set are s⁰ ₁₂ and so on; and filtering eachof these data sets by convolution with a suitable kernel to obtain theelements s_(pq) of the n×n structural matrix S.

[0030] Alternatively step (c) comprises filtering for each point thepartial derivatives u_(xi) by convolution with a suitable kernel toobtain regularized partial derivatives u^(r) _(xi); calculating for eachpoint a symmetric n×n structural precursor matrix S⁰, wherein theelements s⁰ _(pq) equal u^(r) _(xp)·u^(r) _(xq); creating (½)n(n+1) datasets, wherein the elements of the first data set are s⁰ ₁₁ pertaining toeach point, the elements of the second data set are s⁰ ₁₂ and so on; andfiltering each of these data sets by convolution with a suitable kernelto obtain the elements s_(pq) of the n×n structural matrix S.

[0031] The above described alternatives for step (c), wherein aprecursor matrix is calculated allows defining the edge preservationparameter as follows:

[0032] ε=trace(S⁰S)/(trace(S⁰)·trace(S)). In this way the edgepreservation parameter is calculated for each point.

[0033] The present invention also relates to a method of processing animage in order to display the edge preservation parameter. This methodcomprises the steps of

[0034] (a) obtaining an n-dimensional initial image data set, whereineach element u(m0) of the data set is the initial image intensity of apoint of the image, and wherein n=2 or n=3;

[0035] (b) calculating for each point the partial derivatives of u(m) inn directions to obtain n derivative data sets of the partial derivativesu_(xi);

[0036] (c) calculating for each point a symmetric n×n structural matrixS⁰, wherein the elements s⁰ _(pq) equal u_(xp)·u_(xq); creating(½)n(n+1) data sets, wherein the elements of the first data set are s⁰₁₁ pertaining to each point, the elements of the second data set are s⁰₁₂ and so on; and filtering each of these data sets by convolution witha suitable kernel to obtain the elements s_(pq) of the n×n structuralmatrix S; and

[0037] (d) calculating for each point ε=trace(S⁰S)/(trace(S⁰)·trace(S))and attributing the calculated value of ε to each point to obtain animage in which faults are highlighted.

[0038] Alternatively, the method of processing an image comprises thesteps of

[0039] (a) obtaining an n-dimensional initial image data set, whereineach element u(m=0) of the data set is the initial image intensity of apoint of the image, and wherein n=2 or n=3;

[0040] (b) calculating for each point the partial derivatives of u(m) inn directions to obtain n derivative data sets of the partial derivativesu_(xi);

[0041] (c) filtering for each point the partial derivatives u_(xi) byconvolution with a suitable kernel to obtain regularized partialderivatives u^(r) _(xi); calculating for each point a symmetric n×nstructural precursor matrix S⁰, wherein the elements s⁰ _(pq) equalu^(r) _(xp)·u^(r) _(xq); creating (½)n(n+1) data sets, wherein theelements of the first data set are s⁰ ₁₁ pertaining to each point, theelements of the second data set are s⁰ ₁₂ and so on; and filtering eachof these data sets by convolution with a suitable kernel to obtain theelements s_(pq) of the n×n structural matrix S;

[0042] (d) calculating for each point ε=trace(S⁰S)/(trace(S⁰)·trace(S))and attributing the calculated value of ε to each point to obtain animage in which faults are highlighted.

[0043] In the paper ‘Coherence-Enhancing Diffusion Filtering’ by JWeickert, International Journal of Computer Vision vol 31(2/3), pages111-127, (1999), modifications of the method known from theabove-mentioned thesis are discussed. The paper discloses a method forenhancing coherence in images with flow-like structures. This method isnot a method for preserving edges, and no parameter corresponding to theparameter ε of the present invention is disclosed.

[0044] The present invention will now be described by way of example inmore details. At first the equations for a two-dimensional image arediscussed, then for a three dimensional image and thereafter alternativeembodiments.

[0045] In two dimensions, the initial image data set comprises k₁×k₂elements (or points) u(m=0), wherein the value of each element u is theinitial image intensity of a point of the image, and wherein k₁ and k₂are the number of points in the two directions x₁ and x₂ of thetwo-dimensional image. The word pixel is used to refer to such a pointin two dimensions, and in three dimensions, the point is called a voxel.

[0046] The next step comprises calculating the partial derivatives ofall points u(m) in n directions to obtain n derivative data sets of thepartial derivatives u_(xi). The first data set comprises k₁×k₂ values ofu_(x1) (=∂u/∂x₁) and the second data set comprises k₁×k₂ values ofu_(x2) (=∂u/∂x₂).

[0047] Then for each point a symmetric n×n structural matrix S iscalculated, wherein the elements s_(pq) equal u_(xp)·u_(xq). Thestructural matrix S is in two dimensions $\begin{pmatrix}{u_{x1} \cdot u_{x1}} & {u_{x1} \cdot u_{x2}} \\{u_{x2} \cdot u_{x1}} & {u_{x2} \cdot u_{x2}}\end{pmatrix}.$

[0048] Then for each point an n×n diffusion matrix D is calculated,wherein the elements d_(ij) are a function of the elements s_(pq), whichelements s_(pq) are a function of the partial derivatives u_(xp) andu_(xq).

[0049] Then the image comprising the elements u(m=0) is processed in Msteps to obtain a processed image comprising elements u(m=M), whereinafter each step the integer m is increased with 1. The k₁×k₂ valuesu(m+1) of each image are obtained from the values u(m) from a previousimage by solving a diffusion type equation:u(m+1)=u(m)+c.div(ε.Dgrad(u(m))). The number of times the image isprocessed, M, is a predetermined value. In order to remove noise, M issuitably in the range of from 2 to 4, and when in addition thereto theimage has to be simplified, M is suitably in the range of from 5 to 20.

[0050] In this equation, c is a predetermined constant, 0≦c≦1, and ε isa scalar, 0≦ε≦1, wherein ε is close to zero when near an edge and closeto 1 when far away from an edge. The scalar ε can be understood to be afault highlighter.

[0051] In two dimensions,${{div}\left( {ɛ \cdot {{Dgrad}(u)}} \right)} = {\sum\limits_{i = 1}^{2}\quad {\frac{\partial\quad}{\partial x_{i}}{\left( {ɛ\left( {{d_{i1}u_{x1}} + {d_{i2}u_{x2}}} \right)} \right).}}}$

[0052] In three dimensions, the image data set comprises k₁×k₂×k₃ imageintensities u(m). The first derivative data set comprises k₁×k₂×k₃values of u_(x1) (=∂u/∂x₁), the second derivative data set comprisesk₁×k₂×k₃ values of u_(x2) (=∂u/∂x₂), and the third derivative data setcomprises k₁×k₂×k₃ values of u_(x3) (=∂u/∂x₃). The n×n structural matrixS is $\begin{pmatrix}{u_{x1} \cdot u_{x1}} & {u_{x1} \cdot u_{x2}} & {u_{x1} \cdot u_{x3}} \\{u_{x2} \cdot u_{x1}} & {u_{x2} \cdot u_{x2}} & {u_{x2} \cdot u_{x3}} \\{u_{x3} \cdot u_{x1}} & {u_{x3} \cdot u_{x2}} & {u_{x3} \cdot u_{x3}}\end{pmatrix}.$

[0053] In three dimensions,${{div}\left( {ɛ \cdot {{Dgrad}(u)}} \right)} = {\sum\limits_{i = 1}^{3}\quad {\frac{\partial\quad}{\partial x_{i}}{\left( {ɛ\left( {{d_{i1}u_{x1}} + {d_{i2}u_{x2}} + {d_{i3}u_{x3}}} \right)} \right).}}}$

[0054] In the above, the diffusion matrix is not calculated directlyfrom the image intensities. However, suitably the diffusion matrix iscalculated from filtered data. Therefore the step of calculating thestructural matrix suitably comprises calculating for each point asymmetric n×n structural precursor matrix S⁰, wherein the elements s⁰_(pq) equal u_(xp)·u_(xq); creating (½)n(n+1) data sets, wherein theelements of the first data set are s⁰ ₁₁ pertaining to each point, theelements of the second data set are s⁰ ₁₂ and so on; and filtering eachof these data sets by convolution with a suitable kernel to obtain theelements s_(pq) of the n×n structural matrix S.

[0055] In an alternative embodiment, the n×n structural matrix S iscalculated from filtered partial derivatives. To this end, the partialderivatives u_(xi) are filtered for each point by convolution with asuitable kernel to obtain regularized partial derivatives u^(r) _(xi).For each point a symmetric n×n structural precursor matrix S⁰ iscalculated, wherein the elements s⁰ _(pq) equal u^(r) _(xp)·u^(r) _(xq).Then the (½)n(n+1) data sets are created, wherein the elements of thefirst data set are s⁰ ₁₁ pertaining to each point, the elements of thesecond data set are s⁰ ₁₂ and so on. And these data sets are filtered byconvolution with a suitable kernel to obtain the elements s_(pq) of then×n structural matrix S.

[0056] Suitably the edge preservation parameter ε is derived from theimage data set by calculations. Applicant had found that the structuralprecursor matrix S⁰ can be used to estimate the scalar ε, suitablyε=trace(S⁰S)/(trace(S⁰)·trace(S)), wherein trace(A) is the sum of thediagonal elements a_(kk) of the matrix A. In this way ε becomes afunction of the partial derivatives of u(m) in the two or threedirections. And the scalar ε is a fault highlighter that is calculatedfor each point. The scalar is about 1 in the absence of a fault and muchsmaller than 1 in the presence of a fault.

[0057] Suitably step (d) comprises determining the n eigenvalues λ_(i)and n eigenvectors v_(i) of each of the structural matrices S; sortingthe eigenvalues so that λ₁≧λ₂(≧λ₃) and calculating for each point an n×ndiffusion matrix D, wherein the elements d_(pq) equal v_(2p)·v_(2q)(wherein v₂ is the eigenvector pertaining to the smallest eigenvalue ifn=2) or wherein the elements d_(pq) equal v_(2p)·v_(2q)+v_(3p)·v_(3q)(wherein v₂ and v₃ are the eigenvectors pertaining to the smallesteigenvalues if n=3). The symbols v_(ip) and v_(iq) denote the p-th andq-th element of the eigenvector v_(i). As a result the eigenvectorpertaining to the largest eigenvalue does not contribute to thediffusion matrix, and the eigenvector(s) pertaining to the smallereigenvalue(s) do contribute. This inhibits diffusion of the imageluminance in the direction of the eigenvector pertaining to the largesteigenvalue, which latter eigenvector is directed perpendicular to thereflection.

[0058] The suitable kernel is a kernel for low pass filtering, such akernel is symmetric and positive everywhere. A suitable kernel is theGaussian kernel of width σ>0, which is given by the following equation:The Gaussian kernel is used as a${K_{\sigma}(x)} = {\frac{1}{\left( {2\quad \pi \quad \sigma^{2}} \right)^{n/2}} \cdot {\exp \left( \frac{- {x}^{2}}{2\quad \sigma^{2}} \right)}}$

[0059] convolution mask, wherein x is the position of the centre of theconvolution mask and

is the Euclidean norm, and n is the dimension.

[0060] The invention will now be described by way of example withreference to the accompanying drawings, wherein

[0061]FIG. 1 shows the original image;

[0062]FIG. 2 shows the image treated according to the present inventionwith three anisotropic diffusion steps;

[0063]FIG. 3 shows the image treated according to the present inventionwith ten anisotropic diffusion steps;

[0064]FIG. 4 shows the image treated not according to the presentinvention with three anisotropic diffusion steps with no edgepreservation, ε=1;

[0065]FIG. 5 shows the image treated not according to the presentinvention with ten anisotropic diffusion steps with no edgepreservation, ε=1;

[0066]FIG. 6 shows the image treated not according to the invention withisotropic diffusion without edge preservation;

[0067]FIG. 7 shows the image treated not according to the invention withisotropic diffusion with edge preservation according to Perona-Malikwith λ=10; and

[0068]FIG. 8 shows the image of FIG. 1 treated so as to highlight thefaults.

[0069] Reference is made to FIG. 1, which shows the original image. Theimage is 128 pixels by 128 pixels having a value in the range from 0 to255. The picture shows reflectors 1, 2, 3, 4 and 5 and faults 10, 11,12, 14 and 15. For the sake of clarity, not all reflectors and faultshave been referred to with a reference numeral.

[0070]FIGS. 2 and 3 show the original image treated according to thepresent invention. The structural tensor had been calculated as follows,at first for each point the partial derivatives u_(xi) were filtered byconvolution with a Gaussian kernel to obtain regularized partialderivatives u^(r) _(xi); then for each point a symmetric n×n structuralprecursor matrix S⁰ was calculated, wherein the elements s⁰ _(pq) equalu^(r) _(xp)·u^(r) _(xq); creating (½)n(n+1) data sets, wherein theelements of the first data set are s⁰ ₁₁ pertaining to each point, theelements of the second data set are s⁰ ₁₂ and so on. To obtain theelements s_(pq) of the n×n structural matrix S each of these data setsis filtered by convolution with a Gaussian kernel.

[0071] In order to get FIG. 2, three anisotropic diffusion steps wereapplied. The edges are made clearer whilst noise is reduced.

[0072]FIG. 3 shows a smoother picture, with less noise, obtained withthe method according to the present invention after ten anisotropicdiffusion steps.

[0073] In order to show the improvement obtained with the methodaccording to the present invention, reference is made to FIGS. 4-7,which have been obtained with known methods. FIG. 4 shows the imagetreated not according to the present invention with three anisotropicdiffusion steps with no edge preservation, in which ε=1.

[0074]FIG. 5 shows the image treated not according to the presentinvention with ten anisotropic diffusion steps with no edgepreservation, ε=1.

[0075]FIG. 6 shows the image treated not according to the invention withisotropic diffusion without edge preservation; and FIG. 7 shows theimage treated not according to the invention with isotropic diffusionwith edge preservation according to Perona-Malik with λ=10.

[0076]FIG. 8 shows an image.obtained from FIG. 1 processed to highlightthe faults. The method of processing the image comprised calculating foreach point of the original image (FIG. 1) the partial derivatives ofu(m) in n directions to obtain n derivative data sets of the partialderivatives u_(xi); filtering for each point the partial derivativesu_(xi) by convolution with a suitable kernel to obtain regularizedpartial derivatives u^(r) _(xi); calculating for each point a symmetricn×n structural precursor matrix S⁰, wherein the elements s⁰ _(pq) equalu^(r) _(xp)·u^(r) _(xq); creating (½)n(n+1) data sets, wherein theelements of the first data set are s⁰ ₁₁ pertaining to each point, theelements of the second data set are s⁰ ₁₂ and so on; and filtering eachof these data sets by convolution with a suitable kernel to obtain theelements s_(pq) of the n×n structural matrix S. Then for each pointε=trace(S⁰S)/(trace(S⁰)·trace(S)) was calculated, and the value of ε wasattributed to each point to obtain FIG. 8 in which faults arehighlighted in black.

1. Method of processing an image comprising the steps of (a) obtainingan n-dimensional initial image data set, wherein each element u(m=0) ofthe data set is the initial image intensity of a point of the image, andwherein n=2 or n=3; (b) calculating for each point the partialderivatives of u(m) in n directions to obtain n derivative data sets ofthe partial derivatives u_(xi); (c) calculating for each point asymmetric n×n structural matrix S, wherein the elements s_(pq) equalu_(xp)·u_(xq); (d) calculating for each point an n×n diffusion matrix D,wherein the elements d_(ij) are a function of the elements s_(pq); (e)calculating for each point u(m+1) from u(m) using the following equationu(m+1)=u(m)+c.div(ε.Dgrad(u(m))), wherein c is a predetermined constant,0≦c≦1, and wherein ε is a scalar, 0≦ε≦1, wherein ε is close to zero whennear an edge and close to 1 when far away from an edge; and (f)repeating steps (b) through (e) M times to obtain the processed image.2. The method according to claim 1, wherein step (c) comprisescalculating for each point a symmetric n×n structural precursor matrixS⁰, wherein the elements s⁰ _(pq) equal u_(xp)·u_(xq); creating(½)n(n+1) data sets, wherein the elements of the first data set are s⁰₁₁ pertaining to each point, the elements of the second data set are s⁰₁₂ and so on; and filtering each of these data sets by convolution witha suitable kernel to obtain the elements s_(pq) of the n×n structuralmatrix S.
 3. The method according to claim 1, wherein step (c) comprisesfiltering for each point the partial derivatives u_(xi) by convolutionwith a suitable kernel to obtain regularized partial derivatives u^(r)_(xi); calculating for each point a symmetric n×n structural precursormatrix S⁰, wherein the elements s⁰ _(pq) equal u^(r) _(xp)·u^(r) _(xq);creating (½)n(n+1) data sets, wherein the elements of the first data setare s⁰ ₁₁ pertaining to each point, the elements of the second data setare s⁰ ₁₂ and so on; and filtering each of these data sets byconvolution with a suitable kernel to obtain the elements s_(pq) of then×n structural matrix S.
 4. The method according to claim 2 or 3,wherein ε=trace(S⁰S)/(trace(S⁰)·trace(S)).
 5. The method according anyone of the claims 1-4, wherein step (d) comprises determining the neigenvalues λ_(i) and n eigenvectors v_(i) of each of the n×n structuralmatrices S; sorting the eigenvalues so that λ₁≧λ₂(≧λ₃) and calculatingfor each point an n×n diffusion matrix D, wherein the elements d_(pq)equal v_(2p)·v_(2q) (if n=2) or wherein the elements d_(pq) equalv_(2p)·v_(2q)+v_(3p)·v_(3q) (if n=3).
 6. A method of processing an imagecomprising the steps of (a) obtaining an n-dimensional initial imagedata set, wherein each element u(m=0) of the data set is the initialimage intensity of a point of the image, and wherein n=2 or n=3; (b)calculating for each point the partial derivatives of u(m) in ndirections to obtain n derivative data sets of the partial derivativesu_(xi); (c) calculating for each point a symmetric n×n structural matrixS⁰, wherein the elements s⁰ _(pq) equal u_(p)·u_(q); creating (½)n(n+1)data sets, wherein the elements of the first data set are s⁰ ₁₁pertaining to each point, the elements of the second data set are s⁰ ₁₂and so on; and filtering each of these data sets by convolution with asuitable kernel to obtain the elements s_(pq) of the n×n structuralmatrix S; and (d) calculating for each pointε=trace(S⁰S)/(trace(S⁰)·trace(S)) and attributing the calculated valueof ε to each point to obtain an image in which faults are highlighted.7. A method of processing an image comprising the steps of (a) obtainingan n-dimensional initial image data set, wherein each element u(m=0) ofthe data set is the initial image intensity of a point of the image, andwherein n=2 or n=3; (b) calculating for each point the partialderivatives of u(m) in n directions to obtain n derivative data sets ofthe partial derivatives u_(xi); (c) filtering for each point the partialderivatives u_(xi) by convolution with a suitable kernel to obtainregularized partial derivatives u^(r) _(xi); calculating for each pointa symmetric n×n structural precursor matrix S⁰, wherein the elements s⁰_(pq) equal u^(r) _(xp)·u^(r) _(xq); creating (½)n(n+1) data sets,wherein the elements of the first data set are s⁰ ₁₁l pertaining to eachpoint, the elements of the second data set are s⁰ ₁₂ and so on; andfiltering each of these data sets by convolution with a suitable kernelto obtain the elements s_(pq) of the n×n structural matrix S; and (d)calculating for each point c=trace(S⁰S)/(trace(S⁰)·trace(S)) andattributing the calculated value of ε to each point to obtain an imagein which faults are highlighted.